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We study both analytically and numerically the existence, uniqueness, and stability of vortex 
and dipole vector solitons in a saturable nonlinear medium in (2+1) dimensions. We construct 
perturbation series expansions for the vortex and dipole vector solitons near the bifurcation point 
where the vortex and dipole components are small. We show that both solutions uniquely bifurcate 
from the same bifurcation point. We also prove that both vortex and dipole vector solitons are 
linearly stable in the neighborhood of the bifurcation point. Far from the bifurcation point, the family 
of vortex solitons becomes linearly unstable via oscillatory instabilities, while the family of dipole 
solitons remains stable in the entire domain of existence. In addition, we show that an unstable 
vortex soliton breaks up either into a rotating dipole soliton or into two rotating fundamental 
solitons. 



PACS: 42.65.Tg, 05.45.Yv. 



1 



I. INTRODUCTION 



Spatial solitons have been a subject of many studies since their first theoretical prediction [1] . The previous research 

on spatial solitons was driven by their promising applications in all-optical devices in which the light guides and steers 
the light itself [2]. Early works studied optical materials with the Kerr (cubic) nonlinearity which exhibit stable 
fundamental (single-hump) solitons in one spatial dimension [3-5] and collapse in two and three spatial dimensions 
[6,7]. Later works focused on optical materials with saturable nonlinear response such as photorefractive crystals (see 
[8] and references therein). The nonlinearity saturation suppresses the collapse of fundamental solitons in two and 
three dimensions [9-11], which opens the door for their experimental observation in multi-dimensional optical beams. 
The instability of higher-order (multi-hump) solitons is not suppressed by the nonlinearity saturation however [12-16]. 
Internal modes of fundamental solitons in saturable optical materials have also been reported [16,17]. These modes 
are responsible for long-lived shape oscillations. 

Recently, incoherent coupling of spatial solitons in photorefractive crystals was proposed and experimentally demon- 
strated [11,18]. The mutual trapping of incoherent optical beams leads to many novel spatial solitons such as vortex 
and dipole vector solitons [19 28]. The incoherently coupled spatial solitons arc described by a system of coupled 
nonlinear Schrodinger (NLS) equations. A similar system of equations also describes temporal solitons in birefringent 
optical fibers and wavelength-division- multiplexed systems [29-34] . Additionally, vortex vector solitons are known in 
the Bose-Einstein condensation guided by a magnetic trap [35]. 

Vortex and dipole vector solitons in saturable optical materials are interesting for both physical and mathematical 
reasons. Physically, these spatial solitons are novel nonlinear objects. They bifurcate from a coupled state where a 
fundamental soliton in one component guides a small higher-order mode in the other component. Far from the bifur- 
cation threshold, both components strongly trap each other and form a fully coupled vector soliton. Mathematically, 
existence and stability of the vortex and dipole vector solitons are challenging problems due to their complexity. The 
existence of vortex and dipole solitons was established numerically [20,23] and with a heuristic variational method [28]. 
However, analytical expressions for radial-angular dependences of vortex and dipole solitons have not been found. On 
the stability of vortex solitons, the numerical results of [20] suggest that vortex solitons with sm,all vortex components 
are stable and observable for large propagation distances, while vortex solitons with large vortex components are 
unstable. Numerical results of [23] show however that all vortex solitons are linearly unstable, and the instability 
leads to the breakup of vortex solitons into rotating dipole solitons. The discrepancy between [20] and [23] raises 
an open question: are vortex vector solitons with small vortex components really stable or not? On the other hand, 
dipole solitons were found numerically to be always stable in [23]. 

In this paper, we clarify the issues of existence and stability of vortex and dipole vector solitons in a saturable 
nonlinear medium such as photorefractive crystals. First, we address the existence and uniqueness of vortex and 
dipole solitons with the perturbation technique [36,37]. We derive perturbation series expansions for the vortex and 
dipole solitons near the bifurcation point where the vortex and dipole components are small. The analytical formulas 
for these solitons are in excellent agreement with our numerical results. We also prove that those vector solitons are 
unique up to phase-, translation-, and rotation-invarianccs. Next, we study the linear stability of vortex and dipole 
solitons with both the spectral analysis and numerical methods. We show that dipole solitons are linearly stable in the 
entire existence domain, while the vortex solitons with large vortex components are linearly unstable, in agreement 
with [23]. However, we prove that vortex solitons with small vortex components are linearly stable^ confirming the 
results of [20], not [23]. Lastly, we study the nonlinear evolution of linearly unstable vortex solitons. We show that 
an unstable vortex soliton breaks up into a rotating dipole soliton only when the vortex component is below a certain 
threshold. Above this threshold, an unstable vortex soliton breaks up into two fundamental vector solitons instead. 



II. EXISTENCE AND UNIQUENESS OF VORTEX AND DIPOLE SOLITONS 



The mathematical model for two incoherently-coupled laser beams in a photorefractive crystal is well known (see, 
e.g., [20,23]). After variable rescalings, the model can be written as a system of coupled equations: 
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where A is the two-dimensional Laplacian, and s is the saturation parameter. 
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Vector solitons in this model take the form: 



El = u{x,y)e", E2 = w{x,y)e'^', (3) 

where the frequency of the Ei wave has been normalized to one, and the frequency of the £'2 wave is A. The amplitude 
functions u{x, y) and w{x, y) satisfy the nonlinear boundary-value problem with zero boundary conditions on the {x, y) 
plane: 

*!!±M) =„, (4) 

1 + s{\uY + \wY) 
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l + s(|u|2 + |«;|2) ^ ' 

The system (4) (5) may have several types of vector solitons localized in two dimensions. The fundamental vector 
solitons take the form: 

u = c„$(r), w = c^$(r), (6) 

where $(r) is a real-valued, singlc-hump function, r = \/ x'^ + 2/^, and Cu and arc arbitrary complex parameters 
constrained by the relation + |cu,p = 1. These solitons exists at A = 1 where the system (4)-(5) reduces to a 
scalar equation for $(r), see Eq. (10) below. The vortex vector solitons take a general form: 

M = $„(r)e»^ w = <^^{r) e'""^ , (7) 

where n and m are topological charges of vortices in the u and w components, $u(r) and $u,(r) are real- valued 
functions, and (r, 6) are the polar coordinates on the {x, y) plane. The simplest vortex soliton has n = and m = 1 
[20,23]. The multi-pole vector solitons take yet another general form: 

u = U{r,e), w = W{r,6), (8) 

where U (r, 9) and W{r, 9) are real-valued functions and may have multi-hump profiles on the {x, y) plane. The multi- 
pole solitons with a single hump for u{x,y) and multiple humps for w{x,y) were approximated in the variational 
approach [28] by the ansatz: 

u = U{r), w = W{r) cos m9, (9) 

where the number of humps in the w component is 2m. The simplest multi-pole soliton is a dipole soliton which has 
m = 1. 

In this paper, we study the simplest vortex and multi-pole vector solitons. We will refer to them simply as the vortex 

soliton and dipole soliton hereafter. The existence and uniqueness of the vortex and dipole solitons can be studied by 
perturbation methods. The perturbation series expansions are derived in the neighborhood of the bifurcation value 
A = Ao(s), where the w component is small. With the perturbation arguments, we show that the bifurcation value 
Ao(s) is the same for both branches of vortex and dipole solitons, and these solitons bifurcate uniquely from A = Ao(s). 
Setting w = in (4), we find the nonlinear boundary- value problem for the scalar soliton u = Uo{r): 

u'^ + -u'o -U0 + = 0, (10) 

r 1 + suq 

where uo{r) is a real-valued function. We take uo{r) to be the fundamental soliton, i.e., uo{r) > for finite r > 0. If 
w is small in (5), we get a linear eigenvalue problem for the first-order correction w = Wi{x, y): 

Awi - \wx + - — 2_^w;i = 0, (11) 

1 + SMg 

where w\{x,y) is a complex-valued function in general and A is the eigenvalue. The linear equation (11) supports 
localized solutions of the form (/)(r)e^*'"^ at a set of discrete values of A. Since we study the simplest vortex and 
dipole solitons, we set m = 1, and require (/)(r) to be a non-negative function for r > 0. The corresponding eigenvalue 
Ao(s) and eigenfunction (j){r) satisfy the following reduced equation: 
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The eigenvalue A is unique once s is fixed. We normalize the cigcnfunction (A(r) such that it has a maximum value 
one. Numerically, we compute Ao(s) and (j){r) by the shooting method. Fig. la shows the dependence of Aq versus 
saturation parameter s. Fig. lb shows the scalar soliton Uo{r) and the normalized eigenfunction 0(r) at s = 0.5, 
where Aq = 0.2622. 

When (j){r) and Ao(s) are known, a general solution for •Wi{r, 6) can be written in the form: 

wi{r,e) = 4>{r)(coiie + ipii\ne), (13) 

where p is an arbitrary real parameter. In this general solution, we have removed arbitrary rotations and translations 
on the (a;, y) plane as well as an arbitrary phase shift in the w component. We note that the solution (13) is identical 
to the variational ansatz in [28]. 

Below, we use the perturbative method and show that there are only two continuations of the solution (13): for 
p = ±1 and p = Q. When p = ±1, the perturbation series expansion recovers the vortex soliton (7) with n = 
and m = ±1. When p = 0, the perturbation expansion recovers the dipole soliton (8) with a single hump for u 
and a double hump for w. For given values of s and A, the two solutions are unique up to phase-, translation- and 
rotation-invariances. At other values of p, the solution with w^s leading-order term as (13) can not exist. 

The perturbation series expansions for vector solitons in system (4) and (5) take the form: 

u = uo{r) + e\2{r, 9) + eV(r, 6) + 0{e^), (14) 



w = ewi{r, 9) + €^W3{r, 9) + e^w^ir, 9) + O(e^), (15) 

and 

A = Ao(s) + e^\2{s) + e^A4(s) + 0{e^), (16) 

where e is a small parameter, uo{r) is the scalar fundamental soliton solving (10), w\{r, 9) is the first-order correction 
in the form (13), and the cut-off frequency Ao is the eigenvalue of (12). The objective of the perturbation analysis is 
to uniquely determine the coefficients A2, A4, . . . as well as expressions for functions U2, U4, W3, and so on. Once the 
coefficients Aq, A2, . . . have been obtained, we can compute e from A in the expansion (16). Once e is found, together 
with functions uq, U2, wi,W3, . . ., we can approximate the vector soliton by the expansions (14) and (15). Below, we 
will carry out the perturbative calculations to order e^. 

Substituting the perturbation series (14), (15) and (16) into the original equations (4) and (5), at order e^, we get 
an inhomogeneous equation for U2- 



ul{2 + sul) ul _ uo\wi 
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where u is the complex conjugate of u. The linearized operator in the left-hand-side of (17) has a non-empty null 
space spanned by three linearly independent localized eigenfunctions: 



u 



2h = «wo(r), 4ft = "o('') cos 9, 4ft = (0 sin 9. (18) 

These eigenfunctions correspond to the phase and translational invariances of solitons in the scalar u equation. The 

right-hand-side of (17) is orthogonal to M2ft because it is real-valued. It is also orthogonal to Wjft ^^d U2ft because it 
has different angular dependence of 1 and cos 2^ (rather than cos^ and sin^). Therefore, up to phase-, translation-, 
and rotation-shifts, a localized solution to (17) is constructed uniquely in the form: 

U2 = i(l+p')ii2o(r) + i(l-p2)u22(r)cos2^, (19) 



where functions U2o{r) and U22(r) satisfy the equations 
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+ ~U22 



4\ , ul{i + sul) 

U22 + -n—. 2\F"22 



We do not know exact analytical expressions for U2Q{r) and U22{r) but can compute them numerically. 
At order e^, we get the equation for as 



A \ , fx |wi|2 + 2woM2\ 



We denote 



2(1 + sul) 



2(1 + SM^) 



2\2 ' 



and rewrite Eq. (22) in the form: 



(21) 



(22) 



(23) 



AW3 - X0W3 + 



A2-(l+p')/H + -(l-p2)/l2 



1 + SWq 



W3 



A2-(l+p2)/ii- (l-/)/j2 



1 



)COS^ 



)sin^ — -(1 — p )/i2^cos30 — -ip{l —p )h2(j)smW 



(24) 



The homogeneous part in (24) supports two linearly independent localized solutions 0(r) cos and (p{r) sin 0. As 
a result, a localized solution of the non-homogeneous equation (24) exists if and only if the following solvability 
conditions are satisfied: 



1 



rcl>' <i A2 - (1 +p')hi{r) - -{l-p')h2{r) } dr = 0, 



(25) 



-'0 



A2 - (1 + p')/ii(r) + 2 (1 - p')/i2(r) }dr = 0. 



(26) 



We will show below that these solvability conditions define only two perturbation series solutions for vector solitons. 
These solutions correspond to the choice p = ±lorp = 0, which produce vortex and dipole solitons respectively. 



A. Vortex solitons 



If p 7^ 0, we eliminate parameter A2 from the system (25)-(26) and find the solvability condition in the form: 



/•OO 

(l-j9^) / r(j)'^h2{r)dr = 0. 
Jo 



(27) 



The integral in Eq. (27) only depends on the parameter s, not on p. We have checked numerically that this integral 
never vanishes for any s. Thus, the condition (27) is satisfied only when p = ±1. In this case, it follows from (13), (19) 
and (24) that wi = c/)(r)e^*^, U2 ~ U2o{r), and = /(r)e^*^. We can continue the perturbation series expansions 
(14)-(16) to higher orders and find that all U2n {n > 0) corrections are only functions of r, and all W2n+i {n > 0) 
corrections have the form g{r)e^^^ . Thus, the perturbation series solution gives a vortex vector soliton (7) with n = 
and m = ±1. 

When p = ±1, we find from Eq. (25) that the coefficient A2 is 



A2 = A2j,(s) 



2 rcfy'^h^dr 
r(f)'^dr 



(28) 



The functional dependence of A2,; versus s is computed from this formula and plotted in Fig. la (dashed line) 
alongside the cut-off frequency Ao(s). Since the (non- negative) function <l){r) is normalized to have maximum one, the 

perturbation parameter e determines the amplitude (maximum) of the vortex component w with error at the order 
of e^, see (13) and (15). The dependence of e versus A and s can be obtained from the perturbation series (16) with 
error of order e^: 
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Wc compare the analytical formula (29) with numerical results for s — 0.5, where Xq = 0.2622 and ~ 0.1010. A 
dashed line in Fig. 2(a) shows the amplitude of the vortex component w computed from (29). Numerically, vortex 
solitons are computed from the original system (4) (5) by the shooting method. The amplitudes of the u and w 
components are also shown in Fig. 2. In Fig. 2(b), a profile of u{x, y) and w{x, y) components for s = 0.5 and A = 0.4 
is shown. We can see from Fig. 2a that the agreement between the analytical predictions and numerical values on the 
amplitudes of the w-component is good over a wide range of A values. Also, the numerically-obtained amplitude of 
the w-component depends linearly on A, which is in agreement with the perturbation series (14) up to error of order 
e^, since oc (A — Aq). 

B. Dipole solitons 

If p = 0, the condition (26) is satisfied, while the condition (25) gives the correction term A2 in the form: 

- ^''^'^ = • 

The functional dependence of \2d versus s is numerically computed and plotted in Fig. 1(a) (solid line). When A2 is 
given by (30), a localized solution w;3(r. 9) of Eq. (24) exists, and this solution is real- valued. We can further show that 
the perturbation series expansions (14) (16) can be successfully continued to higher orders of e, and a dipole-soliton 
solution (8) can be obtained. This solution is real-valued and has the symmetries 

u{-x,y) =u{x,y), u{x,-y) = u{x,y), w{-x,y) = -w{x,y), w{x,-y) = w{x,y). (31) 

Similar to the vortex soliton case, the perturbation parameter e here gives the amplitude (maximum) of the dipole 
component w with accuracy of 0{e^). The formula for e is still (29), but the A2 value is now given by Eq. (30) 
instead of (28). The comparison between the analytical result (29)-(30) and numerical results for dipole solitons is 
shown in Fig. 3a for s = 0.5. In this case, Aq = 0.2622 and X2d = 0.1174. We see again that the agreement between 
numerical and analytical results is very good over a wide range of A values. In Fig. 3b, the profiles of u{x, y) and 
w{x, y) components of a dipole soliton, computed with numerical iteration methods, are displayed. 

We note that in view of Eqs. (28) and (30), the integral in Eq. (27) is actually 



f 

Jo 



r(l>^h2{r)dr = (2A2d - Aa^) / r(l>^dr. (32) 



Inspection of the A2d(s) and X2v{s) curves in Fig. la immediately confirms that the integral (32) is always positive. 
Thus Eq. (27) holds only when p = ±l. 

To summarize, wc have shown that there are only two vector solitons of the system (4)-- (5) that bifurcate from the 
cut-off frequency A — Ao(s). They are either a vortex soliton (7) or a real- valued dipole soliton (8). The solutions 
are determined in terms of perturbation series expansions up to order e*^. Both solitons are unique up to phase, 
position and rotation invariances. The analytical results are confirmed by numerical calculations. Computations of 
the perturbation series expansions prove the existence and uniqueness of vortex and dipole vector solitons observed 
numerically in [20,23,28]. 

III. LINEAR STABILITY OF VORTEX AND DIPOLE SOLITONS 

In this section, we study the linear stability of vortex and dipole vector solitons by spectral analysis, supplemented 
by numerical computations. Since the linearization operators differ for vortex and dipole solitons, we shall treat the 
two cases separately. 
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A. Vortex solitons 



To study the linear stability of the vortex solitons (7) with n = and m = 1, we linearize the system (l)-(2) with 
the perturbation in the form: 

= e" [$„(r) + u+(r)e-'"^+'^^ + M_(r)e*"^+*^] , (33) 

E2 = e'^'+''^ [^w{r) + w+{r)e-'"^+''' + w_{r)e'"^+^'] . (34) 
The linearization problem can be written in the form: 

iuu+ = -u" u\ + [l + — \ u+ - V u+ - Vuu {u+ + U-) - Vuw {'W+ + W-), (35) 



-iau- = -u"_ - -u'_ + il + ^]u- -V U- - Vuu {u+ + U-) - Vuw {w+ + W-), (36) 



1 / (n-1)^' 

iaw+ = -w" w'^ + I A H ^ — \ w+ -V 11)+ - Vuw (w+ + u-) - Vww {w+ + W-), (37) 



r 



-law- = —w_ 



^w'_ + + ^'^y^ ^ W--VW-- Vuw {u+ + U-) - Vww {w+ + W-), (38) 



where 



(1)2 + (1)2 (1)2 (1) (1) (1)2 
^r I Tr _^_u T/ _ ^u^w ^T ^t, 



1 + , ((1)2+5)2,)' - (l+,($2+$^))2' - (l + ,(ci)2+$2^))2' (1 + ,((^2 + $2^))2 " 

The linearized problem can be formulated in the Hamiltonian form [35]: 

Sh , . 6h 
±iau± = -r^, ±iaw± = (39) 
du± dw± 

where h is the energy quadratic form associated with an eigenvector u = (u+jU-jW+jW-)^ and a linearized self- 
adjoint operator C of the right-hand-sides of the system (35)-(38): 



/>oo 

h={u,Cu)=ia / rdr (|«+|2-|«_|2 + |u;+|2_|^i,_|2). 

Jo 



(40) 



The eigenvalue a is defined by the spectrum of the linearized problem (35)-(38) when u(r) is localized as r ^ oo such 
that the integral (40) makes sense. The eigenvalues could he isolated or embedded into continuous spectrum of the 
system (35)-(38). The vortex soliton is linearly unstable if tlicrc exists an eigenvalue cr for some n such that Re(cr) > 0. 
We note that if (a, n, w+,u_,io+,tw_) is a solution of the linear system (35)-(38), so are ((7, — n, u_, u+, «)_, «)+), 
{—a, —n,U-,u+,'W-,w+) and {—a, n,u+, u_, «)+, u)_). Thus, complex unstable eigenvalues a always come in quartets, 
while real and imaginary eigenvalues a always come in pairs. We also note that eigenmodes {a,n,u+,u-,w+,w-) 
and ((T, — n, u-, u+, ■«)+) give the identical perturbation in (33) and (34), so do (— cr, — n, u_, u+, w_, w+) and 
(—a, n, u+, U-, w+, W-). 

According to the stability theory of solitary waves in the system of coupled NLS equations [38], only eigenvalues 

(T with negative or zero values of the energy quadratic form (40) may bifurcate to the domain Rc(A) > 0, leading to 
instabilities. We shall apply this theory and study the spectrum of the linearization problem (35)-(38) at the cut-off 
frequency A = Ao(s), when $„(r) = Uo{r) and $u,(r) = 0. We will show that the vortex soliton is linearly stable in 
the neighborhood of the cut-off frequency Ao(s) < A < Xds), where Ac(s) < 1 is the instability threshold. In the limit 
A — » Ao(s), the linearization problem decomposes into two linear problems: 



±iau+ = -u'i -'-u'++(l + u+ - - (41) 

r V J 1 + suq {1 + suqY 



and 

±iaw± = -w'^ - -w'+ + 2 )w± - , ° ■2 W±- (42) 
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The first linear problem (41) is the stability problem of a scalar fundamental soliton u = Mo(r) in a saturable medium. 
The linear stability of such solitons has been well established (sec [16] for instance), thus unstable eigenvalues a do not 
exist in (41). The continuous spectrum of the system (41) is located at Rc((t) = and |Im((T)| > 1. The continuous 
spectrum is irrelevant for stability of solitary waves in the system of coupled NLS equations, if no embedded eigenvalues 
with negative energy quadratic form (40) exist [38]. The discrete spectrum of (41) consists of isolated eigenvalues a 
such that Re(tT) = and |Im(cr)| < 1, including the zero eigenvalue at n = and n = ±1 with three eigenhmctions 
(18) and three generalized eigenfunctions. Additional eigenvalues for internal modes exist in (41) for Re(cr) = and 
7^ |Im((T)| < 1. These modes have been determined numerically in [16]. Using those numerical results, we have 
found that the energy quadratic form (40) is positive for all internal modes of the system (41). For instance, only one 
internal mode with n = exists and has positive value of ft for s = 0.5 (which corresponds to a; = —0.5 in [16], see 
Fig. 3). We have also checked that no embedded eigenvalues with |Im((T)| > 1 exist in the problem (41) for s = 0.5. 

The second linear problem (42) is uncoupled for m;+ and w_ . Since the operator in the right-hand-side of Eq. (42) 
is self-adjoint, the spectrum of a is purely imaginary, i.e. Re(cr) = 0. The continuous spectrum of (42) is located at 
|Im((T)| > Aq. Its discrete spectrum consists of isolated eigenvalues with Im(cr) > — Aq for w+ and Im(CT) < Aq for W- 
and can be embedded into continuous spectrum. The discrete spectrum of (42) includes two zero eigenvalues at n = 
with eigenfunctions w± = (j){r), two zero eigenvalues at n = ±2 with eigenfunctions w± = 4>{r), and two non-zero 
eigenvalues a = ±i(l — Ao) at n = ±1 with eigenfunctions 'W± = Uo{r). The zero eigenvalues at n = are induced 
by the phase invariance of the E2 equation, while the zero eigenvalues at n = ±2 are induced by the symmetry of 
the uncoupled problem (42). The eigenvalues a = ±i{l — Aq) with n = ±1 are induced by the arbitrary polarizations 
in the fundamental vector soliton (6) at A = Ao(s). The latter eigenvalues result in negative values of the energy 
quadratic form (40): 

/■OO 

h = -{l- Ao) / ul{r)rdr < 0. (43) 

JO 

We have checked numerically that (42) has no other discrete eigenvalues for s = 0.5. 

Applying the stability theory of solitary waves [38], we count eigenvalues of the problems (41) and (42) that produce 
negative and zero values of the energy quadratic form h. At A ~ Ao(s), only two eigenvalues a = — Aq) give the 
negative energy (43). Several zero eigenvalues give zero energy at n = 0, ±1, ±2. However, zero eigenvalues at n = 
and n = ±1 are preserved at A > Ao(s) due to translation, rotation, and complex-phase symmetries of the system 
(l)-(2). Only two zero eigenvalues of the problem (42) at n = ±2 are not preserved by the symmetry and they can 
move out of zero for A > Ao(s). We shall now consider the shift of these negative-energy and zero-energy eigenvalues 
for A near the cut-off frequency Ao(s). 

We show first that the negative-energy eigenvalues a = ±i(l — Aq) never bifurcate off the imaginary axis for 
A > Ao(s) regardless whether they are embedded or not. Indeed, at any A value, the linearization problem (35)-(38) 
has the exact discrete eigenmode 

u+ = 0, u_ = -$^(r), w+ = $„(r), u;_ = (44) 

for n = 1 and a = i{l — A), and 

u+ = -'^wir), u- =0, w+ = 0, w- = $„(r) (45) 

for n = —1 and a = — A). This result contrasts what happens in the system of coupled NLS equations with 
Kerr nonlinearities, where the negative-energy discrete eigenvalues, which are embedded in the continuous spectrum, 
bifurcate to the complex plane and lead to the instability [36,38]. 

We note that the exact eigenmodes (44) and (45) generate an approximate solution of the system (l)-(2): 

El = $„(r)e« - 7$^(r)e'^+*^^ + 0(7^), E2 = $^(r)e^^+^^^ + ^^u{r)e" + 0(7^), (46) 

where 7 is an arbitrary small parameter. Solution (46) is nothing but the original vortex vector soliton under a 
small rotation in the {Ei^E^) plane. This rotation leaves the original system (1) and (2) invariant and is one of the 
symmetries of the present problem. Thus, although solution (46) appears like internal oscillations of vortex solitons, 
this oscillation does not create any energy radiation and is fundamentally different from internal oscillations discussed 
in [16,17,36]. 

We show next that the zero eigenvalues at n = ±2 move to the imaginary axis (as conjugate pairs) as A > Ao(s) 

and do not create any instability. We will use the perturbation series expansions and will present calculations only 
for the case n = 2 (the case n = — 2 is similar). When A is close to Ao(s), we construct an approximate solution to 
the linearization problem (35)-(38) at n = 2 in the form: 
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u± = eu'i\r)+0{e^), w+ = (f>{r) + e^w^^\r) + 0{e^), W- = €^w^^\r) + 0{e^), cr = e^as + 0(e*), (47) 

where e is the same small parameter as in expansions (14)-(16). Substituting (47) into the system (35)-(38), we find 
an exact solution at order e: ^i^-* = u^-* = U22{r), where U22{r) solves the equation (21). At order e^, we need to 

(2) 

solve the non- homogeneous equation for w^. (r) : 



(2)// , 1 (2)n /> 1\ (2) , (2) /\ • ^ , 2(l){uoU20 + U0U22 + 4''^) 

-V -(^° + ;:^j<' + T^-V = (A2-^<.2)<^ . (48) 



The solvability condition for this equation can be simplified by virtue of Eq. (28), and we find that the eigenvalue 
coefficient 172 is given as 

<.. = 2ii>^. (49) 

Utilizing Eq. (32), we see that 

CT2 = 2i[2A2d(s)-A2„(s)], (50) 

whoso imaginary part is positive from Fig. la. Thus, the energy quadratic form of the bifurcated eigenmode (47) and 
(50) (up to the order e^) is negative: 



/>oo 

h = -e^Im(CT2) / rcjy^dr < 0. 
Jo 



(51) 



The analytical eigenvalue formula (47) and (50) at n = 2 is plotted in Fig. 4 versus A for s = 0.5 (dash-dotted line). 
Numerically, we have determined these eigenvalues for s = 0.5 and various values of A, and the results are plotted in 
Fig. 4 (solid line) as well. When A is close to Aq, the analytical formula agrees well with the numerical values. 

We have shown above that the two zero eigenvalues of the system (42) at n = ±2 move to the imaginary axis when 
A > Ao(s), while the two non-zero negative-energy eigenvalues at n = ±1 remain on the imaginary axis. Thus, we 
conclude that vortex solitons are linearly stable near the cut-off frequency A = Ao(s), i.e. vortex solitons with small 
vortex components arc linearly stable. This result confirms the conclusions of [20] and does not support conclusions 
of [23] where all vortex vector solitons were claimed to be linearly unstable. 

Unstable eigenvalues of vortex solitons may appear far away from the cut-off frequency Ao(s). Indeed, the two 
imaginary eigenvalues u for n = ±2 that bifurcate from zero eigenvalues at A > Ao(s) have negative energy (51). When 
these eigenvalues collide with eigenvalues of positive energy or with continuous spectrum, the oscillatory instability 
may arise [35,38]. We confirm this scenario and compute unstable eigenvalues a of the linear system (35)-(38) with 
the numerical shooting method. The unstable eigenvalues are found exactly at n = ±2 and are shown in Fig. 4 for 
s = 0.5. The unstable eigenvalues appear when A > Ac ~ 0.402, where Ac denotes the frequency for onset of instability. 
These results agree with Fig. 3 of Ref. [23] , where the unstable eigenvalues were found from time-integration of the 
linearized equations under random-noise initial perturbations. Fig. 5 displays our numerical solutions for the vortex 
soliton and the corresponding unstable eigenfunction at s = 0.5, A = 0.5 and n = 2. Thus, for the case s = 0.5, 
unstable eigenvalues exist at A > Ac « 0.402, while vortex vector solitons exist at A > Aq « 0.2622, see Fig. 4. In the 
interval Aq < A < Ac, i.e., 0.2622 < A < 0.402 for s = 0.5, unstable eigenvalues do not exist and the vortex solitons 
are linearly stable. 

We conclude this analysis with two remarks. First, it follows from Fig. 4 for s = 0.5 that the eigenvalues a at 
n = ±2 merge into the continuous spectrum at A « 0.396, while unstable eigenvalues appear at A = Ac « 0.402. Our 
numerical results are inconclusive as to what happens in the narrow interval 0.396 < A < 0.402. This problem is 
left open for future studies. And second, when A is further away from the cut-off frequency Ao(s), the vector vortex 
solution bifurcates into scalar vortex solutions with u = and w = $tt,(r)e'^, see [28]. The scalar vortex soliton has 
additional unstable eigenmodes at \n\ ^ 2 which have smaller growth rates (see [15]). We do not study this bifurcation 
where the family of vector vortex solitons terminates, nor the number of unstable eigenvalues of vector vortex solitons 
near this bifurcation. 



B. Dipole soliton 



To study the linear stability of the dipole solitons (8), we linearize the system (l)-(2) with the perturbation: 
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El = e" [U{x, y) + (ur{x, y) + Ui{x, y)) e"^ + {ur{x, y) - Ui{x, y)) e*^] , (52) 

E2 = e'^' [W{x, y) + {wr{x, y) + Wi{x, y)) e"' + {wr{x, y) - Wi{x, y)) e^'] . (53) 

Here Ur,Ui,Wr and Wi are complex functions and are very small. The linearization problem is then written in the 
form: 

iaUi = —AUr + Ur — {V + 2Vuu) Ur — 2VuwWr, (54) 

iaUr = —Aui +Ui — Vui, (55) 

iaWi = —AWr + XWr — {V + 2Vyjw) Wr — 2VuwUr-i (56) 

iaWj. = —Awi + Xwi — Vwi, (57) 

where 

V = - V = - V — — — V 



l + s(?72 + ;^2)' (1 + S(f72 + ;^2))2 ' (1 + s({72 + ^2))2 ' (1 + s(C72 + VF2))2 " 

The linearized problem can be formulated in the same Hamiltonian form (39) with the energy quadratic form [38] : 



h — ia / {urUi + UiUr + WrWi + WiWr) dxdy. (58) 



At the cut-off frequency A ~ Ao(.s), the same analysis as for the vortex solitons shows existence of a pair of eigenvalues 
a = ±i{l — Xq) with negative values of h and a number of zero eigenvalues with zero values of h. We show again that 
the eigenvalues cr = ±i(l — Aq) with negative energy never bifurcate into complex domain for A > Ao(s). Indeed, for 
any A value, the linearization problem (54)-(57) has the exact solution 

Ur = -W{x,y), Ui = W{x,y), Wr = U{x,y), Wi = U{x,y) (59) 

at cr = i{l — A), and 

Ur = W{x,y), Ui = W{x,y), Wr = -U{x,y), Wi = U{x,y) (60) 

at (7 = — A). 

We study the zero eigenvalues of the system (54)-(57) with perturbation series expansions for A > Ao(s): 

a = e^CTa + 0{e^), Ur = cu^^\r, 9) + 0(e3), ^ = 0{€^), 
Wr = «;(°)(r, 9) + e\'^^\r, 9) + 0{e^), Wi = wf\r, 9) + e^wf\r, 9) + 0{e^). (61) 

Here e is the same small parameter as in expansions (14)~(16), while the functions w'fl{r, 9) are linear combinations 
of the eigenfunctions of the null space of the problem (54)-(57) at A = Ao(.s): 

w^^ = Ci0(r) cos 9 + C2(j){r) sin 9, wf^^ = di4>(r) cos 9 + d24>(r) sin 9, (62) 

where ci, C2, rfi, and d2 are constants. Substituting (61) into the system (54)-(57), we find an exact solution at order 
e: 

u'^r^ = CiU2oir) + (ci cos 29 + C2 sin 29) U22ir). (63) 

where U2o{f) and U22(j') solve the problems (20) and (21). At order e^, four solvabihty conditions are needed for 

solving the non-homogeneous equations for Wr (r, 9) and wl {r, 9). Using Eq. (30), we transform the four solvability 
conditions to the form: 

/•OO /"OO POO pOO 

a2C\ = 0, 0-2^2 = 0, icr2C2 / r<^dr = d2 r(f>^h2dr, ia2di / r((?dr = — ci / r^^ {2hi + h2) dr. (64) 

Jo Jo Jo Jo 

If (72 = 0, then ci = d2 = 0, while C2, di are arbitrary constants. Thus, the zero eigenvalue persists in the system 
(54)-(57) for A > Ao(s) with two eigenfunctions Wr = (/'(r)sin6' and Wi = (j){r)cos9. The two eigenfunctions are 
related to the symmetries of the system (l)-(2) with respect to rotation in 9 and shift of the complex phase. If 
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cr2 9^ however, the system (64) has only the trivial solution: ci = C2 = rfi = ^2 = 0. Therefore, the other two zero 
eigenvalues do not bifurcate to the imaginary axis but simply disappear for A > Ao(s). 

We have analytically proved above that the dipole solitons are linearly stable in the neighborhood of the cut-off 
frequency Ao(s). Moreover, contrary to vortex solitons, there are only two eigenvalues of negative energy for A > Ao(s), 
and they remain on the imaginary axis for all A values [see (59) and (60] . Thus, we conjecture that the dipole solitons 
are linearly stable in the whole domain of their existence. This conjecture is in agreement with the numerical work of 
[23]. We again confirm this result by numerical simulations of the system (l)-(2) linearized around the dipole soliton 
(8). For s = 0.5, we have simulated the linearized system for several values of A between A = 0.3 and A = 0.85. We 
did not find any instability in the linearized system. Since A ~ 0.3 is close to the cut-off frequency Aq = 0.2622 and 
A = 0.85 is close to the end frequency A = 1, we conclude that dipole solitons are indeed linearly stable in the whole 
existence interval. 



IV. NONLINEAR EVOLUTION OF PERTURBED VORTEX SOLITONS 

Here we study the nonlinear evolution of perturbed vortex solitons. The unstable vortex soliton under small 
random-noise perturbations was found in [23] to break up into a rotating dipole vector soliton. We will show below 
that such a breakup scenario holds only when the vortex component of the vortex soliton is below certain threshold. 
Above that threshold, unstable vortex solitons break up into two rotating fundamental vector solitons instead. We 
will also show that the vortex solitons with small vortex components are not only linearly but also nonlinearly stable. 

We consider first the nonlinear evolution of linearly stable vortex solitons. For this purpose, we have simulated the 
system (l)-(2) starting with a linearly-stable vortex soliton under various types of small initial perturbations such 
as random-noise and deterministic ones. We have found that the vortex solitons are also nonlinearly stable for all 
small perturbations. To demonstrate, we select s = 0.5 and A = 0.38, where the vortex soliton has been shown to be 
linearly stable (see Fig. 4). As initial perturbations, we choose 

Si(r,0,O) = (l + a)*„(r,A), E2{r,e,Q) = [l + a)<^^{r,\)e'\ (65) 

where a is a small perturbation parameter, that measures amplification of the vortex soliton by a factor 1 -\- a. The 
simulation result with a = 0.05 is shown in Fig. 6. This figure shows that the perturbed vortex soliton persists the 

nonlinear evolution and exhibits little change of shape even after 300 diffraction lengths. This clearly confirms the 
linear and nonlinear stability of the vortex soliton with s = 0.5 and A — 0.38. Other perturbations to this soliton give 
similar evolution results. 

We study next the nonlinear evolution of linearly unstable vortex solitons. We have shown in Sec. Ill A that these 
solitons possess two unstable cigcnmodes (cr, n, w+, zt_, w+, w_) and (a, — n, m_, U+, w+) with n = 2. The cr versus 
A graph is shown in Fig. 4 for s = 0.5, while unstable eigenfunctions (u+, u_, w_) for s =^ 0.5 and A = 0.5 are 
displayed in Fig. 5. However, we recognize that these two unstable eigenmodes are equivalent in view of Eqs. (33) 
and (34). Thus, any small initial perturbation to the vortex soliton is projected onto this unstable eigenmode which 
grows exponentially, while the rest of the initial perturbation disperses away. For convenience, we choose the initial 
perturbation to be exactly this unstable eigenmode, i.e., 

Ei{r,e,Q) = *„(r) -ha(w+(r)e-2^^ +«_(r)e2*^), (66) 

E2{r,9,0)=e'' [$^(r) + a(«;+(r)e-2'« + u)_(^)e2^«)] , (67) 

where a is a small perturbation parameter. The advantage of this special perturbation is that it shortens the distance 
for the breakup of the vortex soliton and reduces the radiation noise in the nonlinear evolution of the perturbed 
solution. 

We have discovered two breakup scenarios of the unstable vortex soliton with the initial perturbation (66)- (67). 
We confirm that unstable vortex solitons with relatively small vortex components indeed break up into a rotating 
dipole soliton, in agreement with [23]. However, when the vortex component increases above a certain threshold, an 
unstable vortex soliton breaks up into two rotating fundamental vector solitons rather than one dipole soliton. For 
example, when s = 0.5 and a = 0.05, the vortex soliton breaks up into a dipole soliton when 0.402 < A < 0.45, and 
into two fundamental vector solitons when A > 0.45. Indeed, when A = 0.45 (where the vortex component is relatively 
small), the time evolution of the perturbed vortex soliton is plotted in Fig. 7. It is seen that this soliton breaks up 
into a rotating dipole soliton. But when A = 0.5 (where the vortex component is bigger), the time evolution is shown 
in Fig. 8. Here two rotating fundamental vector solitons are formed after the breakup of the unstable vortex soliton. 
We have also found that these breakup scenarios are insensitive to the type of initial perturbation imposed, because 
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we have simulated the evolutions with different values of a in (66) and (67) as well as with other forms of initial 
perturbations such as random noise, but the breakup scenarios do not change. To check the numerical accuracy of 
our simulations, we have used more grid points and wider {x, y) intervals and obtained identical results. Furthermore, 
our results conserve energies of the Ei and E2 components very well. 

Intuitively, it is not difficult to understand the above two breakup scenarios of unstable vortex solitons. When the 
vortex component of the vortex soliton is small, the instability (with n = 2) breaks up the vortex (E2) component 
into two weak humps, while it does not significantly affect the single- hump shape of the fundamental {Ei) component 
since Ei's initial amplitude is much higher. During the subsequent evolution, the two humps of the E2 component 
arc too weak to break the Ei component into two pieces, thus the solution relaxes into a dipole soliton instead of two 
fundamental solitons. However, when the vortex component of the vortex soliton is sufficiently large, the fundamental 
component becomes small (see Fig. 2 and [23]). In this case, instability breaks up both the vortex and fundamental 
components into two pieces, and two fundamental solitons are formed then. 

V. SUMMARY AND DISCUSSION 

To summarize, we have studied both analytically and numerically the existence, uniqueness, and stability of vortex 
and dipole vector solitons in saturable optical materials in (2+1) dimensions. We have shown that the analytical 
expressions for vortex and dipole vector solitons can be constructed with perturbation series expansions near the cut- 
off frequency A = Ao(s). We have also shown that only two vector solitons bifurcate from the same cut-off frequency, 
which are vortex and dipole solitons. Furthermore, we have proved that both vortex and dipole solitons are linearly 
stable when the vortex and dipole components are small. As the vortex and dipole components increase, the family 
of vortex vector solitons becomes linearly unstable, while that of dipole vector solitons remains linearly stable in the 
entire existence domain. We have also shown that unstable vortex solitons break up into a rotating dipole soliton only 
when the vortex component is relatively small. When the vortex component crosses a certain threshold, the vortex 
soliton breaks up into two rotating fundamental vector solitons instead. We expect that our results are significant 
not only for studies of spatial vector solitons in a saturable nonlinear medium but also for studies of Bose-Einstein 
condensation. 

In this paper, we have studied only the simplest vortex and dipole vector solitons which bifurcate from the fun- 
damental u and small w components. One natural question to ask is the existence and stability of other vortex 
and multi-pole vector solitons. The perturbation series expansion method developed in this paper is powerful for a 
systematic study of general vortex and multi-pole vector solitons near their bifurcation points. But this problem lies 
outside the scope of the present article. We note, however, that vortex solitons (7) with \n\ > and |m| > exist, 
and they are expected to be linearly always unstable because each component has non-zero charge and is linearly 
unstable by itself [15]. This expectation is consistent with our preliminary numerical simulations on vortex solitons 
with charges such as n = 1 and m = — 1. 

Recently, three-component vortex and dipole vector solitons in a saturable medium have been investigated [39]. 
The authors found that those solitons are linearly unstable provided their total topological charge is non-zero. In 
view of our results in this paper, that conclusion needs modification. We plan to study that system carefully in the 
near future. 



ACKNOWLEDGMENTS 

The authors appreciate helpful discussions with Yu. Kivshar, Z. Musslimani, B. Sandstede and D. Skryabin. The 
work of D.P. was supported in part by the NSERC grant No. 5-36694 and CFI grant No. 5-26773. The work of J.Y. 
was supported in part by the Air Force Office of Scientific Research under contract F49620-99- 1-0174, and by the 
National Science Foundation under grant DMS-9971712. 



12 





FIG. 1. (a) The cut-off frequency Ao (dashed-dotted) and the correction terms \2v (dashed) and \2d (solid) for vortex and 
dipole solitons as a function of s. (b) The scalar wo(r) soliton and the normalized eigenfunction ^(r) at s = 0.5. 




FIG. 2. (a) Amplitudes of the vortex vector solitons obtained analytically (dashed line) and numerically (solid line) for 
s = 0.5 and various frequencies A. (b) A numerical vortex-soliton solution with s = 0.5 and A = 0.4. 



(a) (b) 




FIG. 3. (a) Amplitudes of the dipole vector solitons obtained analytically (dashed line) and numerically (solid line) for 
s = 0.5 and various frequencies A. (b) A numerical dipole-soliton solution with s = 0.5 and A = 0.5. 
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FIG. 4. Eigenvalues a of vortex solitons versus A at s = 0.5 and n = 2. The cut-off frequency Ao is maxked by Solid line: 
Im(cr); dashed line: Re(cr); dash-dotted line: analytical formula (47) and (50). 




FIG. 5. The vortex solitori (left) and its unstable eigonmodo (right) at s = 0.5, A = 0.5 and n = 2. In the right figure, solid 
lines are the real parts of the eigenfunctions, and dashed lines are the imaginary parts. 



Z=0 Z=300 



FIG. 6. Stable evolution of vortex vector solitons with A < Ac under the perturbation (65) with s = 0.5, A = 0.38, and 
a = 0.05. 
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FIG. 7. Break-up of an unstable vortex soliton into a rotating dipole soliton under the perturbation (66)-(67) with s = 0.5, 
A = 0.45, and a = 0.05. 
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FIG. 8. Break-up of an unstable vortex soliton into two fundamental solitons under the perturbation (66)-(67) with s = 0.5, 
A = 0.5, and a = 0.05. 
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